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Abstract. For a reliable fully-relativistic Korringa-Kohn-Rostoker Green function 
method, an accurate solution of the underlying single-site scattering problem is 
necessary. We present an extensive discussion on numerical solutions of the related 
differential equations by means of standard methods for a direct solution and by means 
of integral equations. Our implementation is tested and exemplarily demonstrated for 
a spherically symmetric treatment of a Coulomb potential and for a Mathieu potential 
to cover the full-potential implementation. For the Coulomb potential we include an 
analytic discussion of the asymptotic behaviour of irregular scattering solutions close 
to the origin (r «C I). 
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1. Introduction 

The Korringa-Kohn-Rostoker Green function method (KKR) in the framework of 
density functional theory [1] has been a powerful tool for electronic structure calculations 
in the past decades. Although it is based on the early work of Korringa [2] and Kohn 
and Rostoker [3] in 1947 and 1954, respectively, an active development of the method is 
present also today [4]. Originally, the KKR method was formulated for non-relativistic 
approximations using muffin-tin spheres for the description of the potentials. Signihcant 
extensions of the original method are given by the relativistic KKR method [5, 6], by the 
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full-potential KKR method [7, 8] and by the combination of both, the fully-relativistic 
full-potential KKR method [9]. These developments allow for precise calculations of 
total energies and forces including relativistic effects. This opens new perspectives for 
many applications [4] such as structure optimisations, equations of states and phonons. 

The KKR method is based on the multiple scattering theory, whereas a periodic 
system is separated into disjoint atomic regions which represent scattering centres. Since 
the Green function is constructed from the regular and irregular single-site scattering 
solutions at each lattice site, an accurate numerical solution of the underlying differential 
equations is crucial. Within the non-relativistic full-potential KKR method, it was 
suggested by Drittler [7] , to solve the underlying single-site scattering problem by means 
of a Lippmann-Schwinger equation, i.e. in terms of integral equations. In general, these 
integral equations can be solved iteratively via Born series. Later it was suggested to 
reformulate the original approach by means of integral equations of Volterra type to 
achieve better convergence properties [10, 11]. Huhne et ah [9] showed that the original 
approach of Drittler can be used as well in the relativistic KKR method, where instead of 
the Schrodinger equation, the Dirac equation has to be solved. One of the disadvantages 
of the original implementation is given by the fact that a certain cut-off radius Rq close 
to the origin has to be introduced in order to avoid problems arising from the treatment 
of the irregular single-site scattering solutions. Recently, it was shown by Zeller [12], 
that this approximation might be inconvenient for materials like NiTi and that the cut¬ 
off radius can be chosen arbitrarily small by using an analytical decoupling scheme and 
a subinterval procedure with Chebyshev interpolations in each subinterval. 

In the present paper we follow a different approach and discuss besides the solution 
via integral equations, the direct solution of the fully relativistic full-potential single-site 
scattering problem by means of standard methods for ordinary differential equations. 
We demonstrate that both the regular as well as the irregular scattering solutions can 
be obtained with high accuracy. To demonstrate this, we first discuss the solution of 
the spherical Coulomb potential. Furthermore, we derive the asymptotic behaviour of 
the relativistic irregular scattering solutions and show that within the limit towards 
a non-relativistic description their asymptotic behaviour differs significantly from the 
non-relativistic irregular scattering solutions. Second, we discuss the solution of the 
fully relativistic full-potential single-site scattering problem for a Mathieu potential in 
simple cubic lattice and compare different applied solvers. It will be shown, that linear 
multi-step methods are a reasonable choice, and therefore, Adams-Bashforth-Moulton 
predictor-corrector methods will be used to calculate the band structure by means of 
the Bloch spectral function. An analytically calculated band structure for the Mathieu 
potential and the numerically obtained results will be compared. 
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2. Single-site scattering 


In this section we focus on the Dirac or Kohn-Sham-Dirac equation [13, 14] in atomic 
Rydberg units (h=l, m=|,c=|~ 274) applied to non-magnetic systems, 


• V + 


0A(r) = WMf^, 


( 1 ) 


where the index A denotes a combined index According to Huhne et ah [9] 

the equation above can be solved via expanding the solution 0A(r'), which is a four 
component spinor-function, into so called spin-angular functions 



^A'A(r)XA'(^) 

ifA'A{r)XA'{r) 


( 2 ) 


The spin-angular functions, XA^r) = are given by a linear 

combination of complex spherical harmonics Yl^{r) and the Pauli spinors ^rusi 
whereas the expansion coefficients C*™® denote the Clebsch-Gordan coefficients 
— nis) mg') [15]. Using the orthogonality of the spin-angular functions and by 
defining the matrix V_{r) with components Vaa' = (xa| Pefr |xa') the following system of 
coupled first-order ordinary differential equations can be obtained: 
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By definition, the potential Ves{r) is constructed such that it is non-zero within the 
Wigner-Seitz cell of a certain atomic site and zero outside. To obtain this behaviour 
Ueff(u) is given as the product of a non-spherical potential Uesif) and the shape- 
truncation function 0(f^ [16, 17], 


Uefr(A) = Ueff(f)0(r4. 


(5) 


The elements of the potential matrix of the effective potential 14fr(A), based on an 
expansion in terms of complex spherical harmonics, are given by 
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For the construction of the Green function, it is necessary to obtain two linearly 
independent solutions 4 a(A). First, a solution Ra(A), which is regular at the origin, 
and second, a solution which is singular at the origin. Outside of the Wigner- 

Seitz cell both solutions can be constructed from the spherical Bessel functions ji{r) 
and the spherical Hankel functions hi{r) = ji{r) +ini{r), which are a linear combination 
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of the spherical Bessel and the spherical Neumann functions ni{r). By introducing the 
single-site scattering t-matrix t, the regular single-site scattering solution outside of the 
Wigner-Seitz cell can be written in the following way [8, 9], 


Doutside 

/Ia 


(^1 = 

A' 


ji'ikr)xA'i^ 

eA'jj'ikr)xA'if^ 


6a'a - ik 
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tA'A , 


(7) 


by using the abbreviations ca = ;^|^sign(fi:), I = I — sign(K), and the relativistic 
momentum k = with e = W — c^. The associated irregular single-site 

scattering solution is given by 
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outside 
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( 8 ) 


In the following we will discuss two different approaches for the numerical solution of 
equation (3) and (4) inside the Wigner-Seitz cell. 


2.1. Numerical solution of the differential equations 


For a numerical treatment of equation (3) and (4) it is common to transform the large 
and the small component of the solution according to 

fA’A{r) = and 5 'a'a(’^) = (9) 

cr r 

By introducing the matrices Uf'{r) = ]4(r) — £ and U_~ (r) = l) L^ — ^ YJn) as 
well as the spin-orbit coupling matrix iF = {k'^aa'}; ^^le following differential equations 
in matrix form are obtained, 

I: :^Q(r)= -K-Q{r) + U+{r)-P{r), (10) 

clr= r — = — — 

1 ^ 

II: ^P{r) = --K-P{r) + U-{r)-Q{r). (11) 

dr— r — — — = 

Employing appropriate initial conditions the coupled equations (10) and (11) can be 
solved. After this step, the regular scattering solution has to be normalized at the 
boundary of the Wigner-Seitz cell according to equation (7). The normalized expansion 
coefficients of the regular single-site scattering solution inside the Wigner-Seitz cell 
Pff;^^‘^{r) can be found as a linear combination of the numerically obtained solutions 
PA'Air), 

e(r) = ^ PA'A"(r) aA.A, and gk"l""(r) = Qa'aA^) oa'^a • (12) 

A" A" 


By introducing the matrices X(rBs) = {^"bsf(^’"bs)^aa'} and x(rBs) = 

{rBsceAa:j(fcrBs)^AA'} {x = j, h) at the radius of the circumscribing sphere tbs of the 
Wigner-Seitz cell, the following algebraic system of linear equations can be formulated. 


£(rBs) ikEfrjis) A / ^ ^ 1 {tbs) A 
£(rBs) ikh{rBs) ) \t ) ~ \ |(’^bs) j ' 
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To discuss the numerical solution of the ordinary differential equations (10) and (11) 
we used Matlab [18] and compared the performance of seven different methods, which 
will be introduced briefly in the following. First of all, we used the standard solvers 
odel 13 and odel5s, which are an Adams-Bashforth-Moulton predictor-corrector method 
of variable step size and variable order 1,..., 13 and an implicit method of variable 
step size and variable order 1,...,5 based on the numerical differentiation formulas 
[19, 20]. Furthermore, the Dormand-Prince method [21] [ode^S) and the Bogacki- 
Shampine method [22] {ode23) are used. Both solver are explicit methods of Runge- 
Kutta type with variable step size and orders 5 and 3, respectively. Since the Bulirsch- 
Stoer algorithm is used to solve the single-site scattering problem within some KKR 
codes, we included a Gragg-Bulirsch-Stoer algorithm (GBS) with both adaptive order 
and step size [23]. Last but not least, the fixed step size Runge-Kutta method of order 
4 (RK4) and the fixed step size Adams-Bashforth-Moulton predictor corrector method 
of order 5 (ABB) are taken into account, since both are widely applied within the KKR 
community [24]. For the ABB method, we use an explicit Adams-Bashforth method 
of order 5 as predictor which is corrected by an implicit Adams-Moulton method of 
order 5. The corrector step is repeated n times, until the change of the solution due to 
the corrector step is smaller than a requested accuracy (PE(CE)"' method). For RK4 
and ABB the ordinary differential equation was transformed analytically by assuming a 
logarithmic scale x = log(r). 

2.2. Numerical solution via integral equations 

Besides a direct solution of the differential equations (10) and (11) it is possible to solve 
the single-site scattering problem by means of integral equations. Suppose that the 
single-site scattering solution of the Dirac equation is known for an arbitrary reference 
system with the potential Ve^if). These solutions are denoted by 4/A(r), which could 
be either regular or irregular single-site scattering solutions. By using the Lippmann- 
Schwinger equation [25], 

$^(F) = TA(r) + f d^r f, F) AV{P) <Fa(F) , (14) 

JQws 

it is possible to calculate the solutions <Fa(F) corresponding to the potential V)jfr(r) = 
AV{r) + 14fr(F). The quantity Gg(E;r, F) in equation (14) denotes the single-site 
scattering Green function of the reference system and is constructed from the regular 
as well as the irregular single-site scattering wave functions R;^{E,r) and H\{E,E) in 
the following way: 

jg(E;F,F) = (^flAiE,r)nAiE,Pre{r-r') + ME,r)d^{Eyreir-r')') . 

A 

The quantities R\{E,r\i)^ and H(l{E,rn)^ denote the so-called left-hand side solutions, 
which were discussed in detail by Tamura [26]. By expanding the single-site scattering 
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wave function in terms of the spin-angular functions (2) and by using equation (15), the 
Lippmann-Schwinger equation (14) for the regular and irregular single-site scattering 
wave functions can be reformulated in terms of the following matrix equations, 


RA'K{z]r) = 


A" 


RA'A"{z;r) AA'',A{z;r) + HA'A"{z;r) BA",A{z]r) 


(16) 


and 


HA'A{z;r) = ^ 


A" 


flA'A"(2;r)CA»,A(2;r) + //a.a"( 2; r)-D a",a(2; >■) 


(17) 


The matrices A{z,r), ^{z,r), Cl{z,r) and D.{z,r) contain the integrals over the radial 
amplitudes and as well as the functions and which 

include the solutions of the reference system and the change of the potential AV{r), 


AA',A{^;r) = 6 a',a (18) 



DA',A{^;r) = 6 a',a 



The explicit expressions of and are given by 

= -!p^gJ,j,(2;r')AVX+A,( 
Ai 



(23) 

(24) 

(25) 

(26) 


In general, the Lippmann-Schwinger equation (14) is solved iteratively in terms of Born 
series by starting with the reference solution <h®(r) = TA(r) on the right hand side of 
equation (14) and by calculating a better approximation $^(1^) by integration. This 
scheme is repeated, until the difference between two successive approximations, 
and is reasonably small. 
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3. The Coulomb potential 


3.1. Asymptotic behaviour 

The Coulomb potential is spherically symmetric and hence, the expansion into spherical 
harmonics only consists of one (spherically symmetric) component Voo(r), 

Vimi'^) = - ^—hiphm,o- (27) 

Wtt r 


The associated ordinary differential equations are independent of the relativistic 
magnetic quantum number fv and diagonal in k, but coupled with respect to g^. and 
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(28) 

(29) 


It was pointed out by Swainson and Drake [27] that equations (28) and (29) can be 
decoupled by transforming the radial solutions according to 


Mr) \ ^ ( 1 27 \ / g^{r) A 

Lir) J [X 1 J [ Mr) J 


(30) 


Here, the newly introduced quantities are given by X = and 7 = y By 

using the abbreviations e=^,a;^ = (| — ^)^, and by combining equation (30) together 
with (28) and (29), two differential equations of second order can be obtained, which 
are similar to the form of the radial Schrodinger equation. 
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/«(>■) = 0 . 


(31) 

(32) 


Close to the origin (r 1 ), in the asymptotic limit, only the angular momentum terms 
(~ r“^) have to be taken into account, whereas the potential (~ r“^) and the constant 
factor can be neglected. The solution of the resulting differential equations are 
rational functions of the following form. 


gn = cM + C2r ^ \ (33) 

/k = + C 4 r"^. (34) 

Since both solutions g^ and /«, are linear combinations of g^ and /«,, it can be verified, 
that the leading term for the irregular solutions is given by r~'^~^. In the non-relativistic 
limit (^^ ~ 0 ) the asymptotic behaviour of an irregular s-wave function (k = — 1 , 






















Figure 1: Double-logarithmic plot of the real part of relativistic irregular single-site 
scattering solutions for a Coulomb potential {Z = 79) close to the origin. 

I = 0) close to the origin is given by which is in contradiction to the 

non-relativistic solution = r~^. Since the relativistic quantum number k takes on 
values either k = Iotk = —1 — 1, this disagreement can be generalized for all irregular 
wave functions with k = —I — 1. To verify the solution behaviour for r 1, a double 
logarithmic plot of the numerical solution for Z = 79 and c = ^ and various values for 
K is illustrated in hgure 1. The predicted asymptotic behaviour is clearly revealed. 

3.2. Numerical accuracy 

For the discussion of the numerical accuracy for the solution of the differential equation 
(equations (31) and (32)), a Coulomb potential with an atomic number of Z = 79 
was used since it represents the element gold (An). Due to the large atomic mass 
and non-negligible spin-orbit coupling, which causes the typical golden colour, it is a 
prominent example for relativistic effects. Solutions were obtained up to a maximal 
angular momentum quantum number of / = 5. To obtain a scattering state, the energy 
of the scattering solution was chosen to be £ = 1. In the following test, a spherical 
cell was constructed including a minimal radius of tq = 10“^ and a maximal radius of 
^BS = 3. 

The numerical solutions obtained by using different solvers were compared with a 
reference solution, which was obtained by using odd 13 and very high accuracy goals. 
For the solvers taken from the Matlab ode-suite absolute and relative tolerances were 
chosen to be equal and values between 10“^ and 10“^° were used. The maximal relative 
error of the numerically obtained solutions as a function of the number of right-hand side 
evaluations of the differential equation is shown in hgure 2. First of all, it can be verihed 
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Figure 2: Maximal relative error versus number of right-hand side (RHS) evaluations of 
the Dirac equation for a Coulomb potential using different methods for the solution. 


that the numerical accuracy of the solution of the regular single-site scattering solution 
(hgure 3(a)) is similar to the results of the irregular single-site scattering solution (hgure 
3(b)). Hence, the following statements can be given for both types of solutions. As 
summarized by Soderlind et ah [28] an ordinary differential equation is called stiff, if 
the numerical solution via implicit methods performs much better than the solution 
via explicit methods. However, for the present example it can be verihed that the 
method odellS, which is a method for non-stiff equations, performs better than the 
implicit method odelSs. Since we observed poor performance for other implicit methods, 
we conclude that the underlying differential equations for the Coulomb potential are 
non-stiff. The performance of the Dormand-Prince method {ode^S) is very good for 
crude tolerances and becomes comparable to odelSs for hne tolerances. The Adams- 
Bashforth-Moulton predictor-corrector method of order 5 with hxed step size {AB5) is 
rather expensive for high tolerances. But, due to the higher order, the performance is 
better in comparison to the methods of the Runge-Kutta type ode23 and RK4, if high 
accuracy goals are demanded. Also in comparison to ode45, which has the same order, 
less evaluations of the right-hand side of the differential equation for high accuracy 
goals are necessary; it is therefore more efficient. The implemented Gragg-Bulirsch- 
Stoer method with both adaptive order and step-size (GBS) [23] is able to solve the 
differential equations with very high orders. But, for the accuracy goals in practice 
(«i 10“®) it needs about hve times as many evaluations of the right hand side of the 
differential equation in comparison to odellS. 

In many implementations of the KKR method, solvers with hxed step size are 
used for spherically symmetric atomic potentials, whereas a logarithmic mesh of type 
X = log(r) or similar is employed [24]. The reason for this is that the wave functions 
are highly oscillating close to the nucleus, and are smooth for larger values of r. To 
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(a) Solution on a radial mesh 

Figure 3: Step size for the solution of the 
methods. 



X = ln(r) 

(b) Solution using a logarithmic scale 
Coulomb-Dirac problem using adaptive 


verify that a logarithmic mesh is a reasonable choice, the step size for adaptive methods 
close the origin was investigated (see figures 3(a) and 3(b)). The step size during the 
numerical solution of the Coulomb-Dirac problem using odellS, odel5s, and ode^B is 
illustrated in figure 3(a) and compared with the step size of a logarithmic mesh. It 
can be verified that the step size used by ode^B is similar to the characteristic of a 
logarithmic mesh. However, the stair-case-like behaviour of odellS and odelBs occurs 
due to the step-size strategy of the method itself, i.e. the change of the step-size is 
avoided as much as possible [19]. Analogously, it is possible to transform the differential 
equations (31) and (32) to a logarithmic scale x = log(r) and to investigate the varying 
step size of the methods odellS, odelBs and ode^B for the numerical solution of the 
transformed equations (see figure 3(b)). It can be verified, that all three solvers adopt 
a constant step size close to the origin {x < —5), which again reassures the choice of a 
logarithmic scale for methods with fixed step size. 

4. The Mathieu potential 

AT Evaluation of the potential 

A suitable test system for our full-potential implementation is given by a Mathieu 
potential [29], 

V{f) = —Uq {cos{ax) + cos{ay) + cos{az)) , (35) 

that is periodic and represents here a simple cubic crystal structure. Moreover, it is 
highly anisotropic with respect to different crystallographic directions and thus, it is 
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ideal to serve as a test system. For our method it is necessary to express equations (35) 
in terms of complex spherical harmonics (see section 2). By using the exponential form 
of the cosine function and by applying Bauer’s theorem, 

OO I 

e‘" = t'n(kr)Yr(?)Yr(kr . (36) 

/=0 m=—l 


we end up with the following equation, 

OO I 


V{f^ = -V^Uo E Of) ^2^ + fim] ji{kr)Yr{^, (37) 


/=0 m=—l 

with coefficients fim given by 


flm 


^ f ((_!)? +,„i. 


even 


m is odd 


(38) 


4-2. Numerical Accuracy 

Analogously to section 3.2, where the numerical accuracy was discussed for the Coulomb 
potential, our numerical test environment was used to solve the differential equations 
(10) and (11) for the Mathieu potential. According to Yeh et ah [29], the lattice constant 
was chosen to be a = 27r and the pre-factor Uq was set to Uq = —0.5. The expansion of 
the solution in terms of spin-angular functions was evaluated up to a maximal angular 
momentum quantum number of /max = 5. Due to the high value of /max the matrices P 
and Q in (3) and (4) are of dimension 72 x 72. 

The maximal relative error between a very precise reference solution and the 
numerically obtained solution of different solvers versus the number of evaluations of the 
right-hand side of the differential equations are shown in hgure 4. The general behaviour 
using different methods is similar for regular and irregular single-site scattering solutions. 
It can be verihed that for a particular method and for the same number of evaluations of 
the right-hand side of the differential equation the maximal relative error for the regular 
solution is by approximately 2 orders of magnitude smaller compared to the irregular 
solution. This is due to very large absolute values of the irregular single-site scattering 
solutions close to the origin. As shown in the example for the Coulomb potential, we 
conclude that the underlying differential equations for the Mathieu potential can be 
regarded as non-stiff, since the performance of the method odd 13 is much better than 
the performance of the method odd5s [28]. Again, the performance of the Adams- 
Bashforth-Moulton predictor-corrector method of order 5 {AB5) is worse than the 
explicit Runge-Kutta methods RK4 and ode45 for crude tolerances. However due to 
the higher order it is more efficient if high accuracy goals are required. Especially for 
the irregular solutions ode23 fails to give a reasonable solution for an appropriate amount 
of evaluations of the right-hand side of the differential equations. 
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Figure 4: Maximal relative error versus number of right-hand side (RHS) evaluations of 
the Dirac equation for a Mathieu potential using different methods for the solution. 


The differential equations (10) and (11) are characterized by the effective potential 
and since the Mathieu potential is analytically known, methods with adaptive step- 
size like odellS are a reasonable choice. In general, the effective potential within each 
iteration of the KKR method is given on a discrete mesh and, hence, a method with hxed 
step size is appropriate, since interpolations of the potential between mesh points are 
avoided. For the Mathieu potential, it can be verihed (see hgure 4) that the numerical 
solution of the full-potential single-site scattering problem can be obtained by linear 
multi-step methods for non-stiff equations, e.g. by applying an Adams-Bashforth- 
Moulton predictor-corrector method. Therefore, the method AB5 was implemented 
within the computer code Hutsepot [30]. 

Besides the direct solution of the differential equation, the solution of the Dirac 
equation for a Mathieu potential was investigated by means of integral equations. For 
the integration, the trapezoidal rule as well as the Simpson rule was used. In hgure 5 the 
maximal relative error of the solution is plotted versus the number of iterations within 
the Born series for the solution of the Lippmann-Schwinger equation for an irregular 
single-site scattering solution. It can be verihed that the maximal relative error saturates 
quickly and the Born series converges after three iterations. This is in perfect agreement 
with the observation of Drittler [7] for the non-relativistic method. Since the deviation 
from the exact solution is dominated by the error of the numerical quadrature, the 
accuracy of the solution can be improved by either increasing the number of mesh points 
for the quadrature or by improving the integration scheme, e.g. using the Simpson rule 
instead of the trapezoidal rule. In this way the maximal relative error can be decreased 
by about one order of magnitude for the same number of points, where the integrand is 
evaluated. 
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Figure 5: Maximal relative error versus number of iterations within the Born series 
for the solution of the Dirac equation for an Mathieu potential by means of integral 
equations. 


4-3. Band structure 


After the solution of the Dirac equation at each atomic site n, the regular and 
the irregular single-site scattering solutions R'^{E,rn) and H^{E,rn) can be used to 
construct the relativistic multiple-scattering Green function, 


G{E; r; + Rn, + RJ = \ Rl{E, r-4)g^^ {E)Rt,{E, 

— ^^ =AA' 

A,A' 


(39) 


Analogously to the single-site scattering Green function (15), R^{E,fn)^ and 
H4t{E,fn)^ denote the left-hand side solutions [26]. The matrix (E) is called 

=AA' 

structural Green function matrix and can be obtained from the structure constants 
G^{E) of the free space and the single-site scattering t-matrices T{E) = |f”(F^)} at 
each site n, 

g{E) = g{E) [L - UE) g{E)\ . (40) 

Since the Green function has a singularity if the energy hits an eigenvalue of the 
Dirac Hamiltonian, the band structure of a system can be calculated by means of the 
Bloch spectral function, which is given by the imaginary part of the fc-resolved multiple¬ 
scattering Green function [31], 


41b (fc, W") 


-Im Tr 

TT 



G + Rj I W) 


(41) 


The summation in the last expression belongs to all lattice vectors Rj of the lattice 
denoted by L. The trace is a summation of the diagonal components of the 4x4 
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Figure 6: Comparison of the relativistic Bloch spectral function calculated with Hutsepot 
(black) and a band structure calculated by means of a plane-wave approach (blue) for 
the Mathieu potential. 


dimensional Green function. The great advantage of the Mathieu potential is given 
by the fact, that the band structure can be calculated analytically and, therefore, it is 
a good test system for our numerical implementation. In general, a large number of 
terms are necessary within the angular momentum expansion of the Green function to 
obtain the correct energy bands. For illustrative purposes, the Bloch spectral functions 
obtained for different maximal angular momentum (/max = 2, 3,4) within the expansion 
(39) together with the analytically obtained result are compared in Fig. 6. As can 
be seen, high values of /max are necessary to achieve a reasonable agreement between 
the band structures. However, the numerical result for /max = 4 reflects the general 
behaviour of the analytical band structure quite nicely. 
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5. Conclusion 

An elaborate discussion of the numerical accuracy for the solution of the spherical 
Coulomb potential and the Mathieu potential is presented. The solutions were obtained, 
hrst, by various standard methods for the solution of ordinary differential equations 
and, second, by means of an iterative solution of the associated Lippmann-Schwinger 
equation. From the performance of the solvers for stiff and non-stiff problems we 
conclude that the differential equations are non-stiff for both systems. The numerical 
solution can be obtained by using linear multi-step methods like the Adams-Bashforth- 
Moulton predictor corrector method. For the Coulomb potential the asymptotic 
behaviour close to the origin (r 1) was investigated and it could be shown, that 
the non-relativistic limit of the irregular single-site scattering solution shows a different 
asymptotic behaviour in comparison to the associated non-relativistic irregular single¬ 
site scattering solution for the case k = —I — 1. 
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